Precursor-mediated crystallization process in suspensions of hard spheres 
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We report on a large scale computer simulation study of crystal nucleation in hard spheres. 
Through a combined analysis of real and reciprocal space data, a picture of a two-step crystallization 
process is supported: First dense, amorphous clusters form which then act as precursors for the 
nucleation of well-ordered crystallites. This kind of crystallization process has been previously 
observed in systems that interact via potentials that have an attractive as well as a repulsive part, 
most prominently in protein solutions. In this context the effect has been attributed to the presence 
of metastable fluid-fluid demixing. Our simulations, however, show that a purely repulsive system 
(that has no metastable fluid-fluid coexistence) crystallizes via the same mechanism. 
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The crystallization process in complex fluids is not triv- 
ial. For systems such as solutions of proteins, alkanes, 
and colloids it has been shown that crystal nucleation 
rates can be enhanced considerably if the supersaturated 
liquid is quenched to a state that lies close to a metastable 
fluid- fluid critical point [H-Q- The enhanced nucleation 
rate is generally attributed to the fact that the den- 
sity fluctuations occuring in the vicinity of a metastable 
fluid-fluid critical point enable the system to evolve via 
a two-step process. First dense, amorphous precursors 
form and then the crystallization process takes place in- 
side these. The prerequisite of this process scenario, the 
metastable fluid-fluid critical point, is easily realized in 
the systems listed above, which exhibit an interplay of 
repulsive and attractive interactions. 

However, it is worthwhile asking whether the two-step 
process occurs more generally. Surprisingly, there have 
been experiments indicating two-step crystallization oc- 
curing also in hard sphere systems, the simplest model 
system for liquids and crystals (see e.g. Ref. [9(). As the 
interaction energy between two hard spheres is either zero 
(no overlap) or infinite (overlap) , the phase behaviour of 
the system is purely determined by entropy. In particu- 
lar, for one component hard spheres there exists a stable 
crystalline phase but no metastable fluid-fluid demixing 
region. 

The crystallization kinetics in colloidal hard sphere 
systems has been studied experimentally using predom- 
inantly time resolved light scattering |9l- ll5l| a nd to a 
lesser extend real-space imaging techniques [16 [4l9f . In 
the scattering experiments described in Refs. |9l. HH |20| 
the time-evolution of the structure factor has been inter- 
preted using a two-step process model: In the induction 
stage precursors (compressed, structurally heterogeneous 
clusters) slowly grow. Then the precursors are converted 
into highly ordered crystals in a fast, activated process. 
In Ref. !9] it was suggested that size polydispersity lim- 
ited growth is responsible for the induction stage. How- 



ever, later it was argued that the precursor stage be- 
haves in a similar fashion, regardless of polydispersity or 
of metastability suggesting that the precursor nucleation 
and the following conversion is not a special feature of 
polydisperse samples, but that it might constitute a fun- 
damental process of crystal nucleation [l5j . 

The complete mechanism of precursor to crystal con- 
version is still unknown and difficult to obtain via struc- 
ture factor analysis alone. A real space experiment would 
be highly desired. But as size polydispersity cannot be 
avoided in a lab experiment a high precision computer 
study appears to be the best choice. 

Hard spheres are easily realized on the computer. Sim- 
ulations of hard s phe re crystallization have been reported 
e.g. in Refs. |2ll-l2a| and of crystal nucleation kinetics 
e.g. in Refs. [24l. |25|. One should bear in mind, however, 
that nucleation is a typical example of a rare event, i.e. an 
event that has a low reaction rate but high impact on the 
properties of a system. Rare events, and in particular 
nucleation, are very often simulated by methods that are 
based on transition state theory [2|| . The basic assump- 
tion underlying transition state theory is that the rare 
process can be reduced to the dynamics of "slow" vari- 
ables which evolve in an effective free energy landscape 
formed by the "fast" variables (i.e. those that quickly 
adopt a Boltzmann distribution). The choice of slow 
variables already presupposes a certain dynamic scenario 
which does not necessarily hold in the experimental sys- 
tem. For instance, the transition rates computed in [25[ 
using transition state theory correspond to a one-step 
crystallization process and the results are in disagree- 
ment with the experimental results. Therefore we have 
carried out a computer simulation study which follows 
the nucleation kinetics directly (without using any bias- 
ing scheme that would require underlying assumptions on 
the nucleation pathway) and therefore allows for direct 
comparison to experiments. 

We have performed Monte Carlo simulations of N = 
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216,000 hard spheres in a box of volume V — 59.2 x 
59.8 x 59.2 D 3 , where D is the particle diameter. In the 
following we use the particle diameter as unit of length, 
ksT as unit of energy and attempted MC moves per par- 
ticle ("sweeps") as unit of time. The system was pre- 
pared by a fast pressure quench from the stable liquid. 
The subsequent crystallization dynamics were simulated 
at fixed N, V and T by small translational MC moves 
only - a method which mimicks Brownian dynamics on 
long time-scales [13] ■ We used this approach, because 
it requires relatively little CPU time per move, a neces- 
sary property when running a simulation of such a large 
size. We let the system evolve for 10 6 MC sweeps and 
sampled observables every 5,000 sweeps. The number 
density after the quench was N/V = 1.03 (volume frac- 
tion <p = 0.54), a value which lies in the liquid-solid co- 
existence region close to the density of the solid at co- 
existence (ca. 9% above the coexistence density of the 
liquid). This corresponds to a chemical potential dif- 
ference between the metastable liquid and the stable, 
almost completely crystalline state of A// ~ —0.58 per 
particle. Given that the interfacial tension is of the or- 
der of 0.5 [H, [2^| one would expect the system to be far 
beyond the classical nucleation regime and to crystallize 
almost instantaneously. The self diffusion constant was 

D s = 2.3 ■ 10~ 5 #M Cswecp S - 

We monitored crystallization by means of the q6q6- 
bond order parameter [3(| HH . For a particle i with n(i) 
neighbours, the local orientational structure is character- 
ized by 

n(i) 

qim (i) := -7TT Y] Yim (fij ) , 
»(*) 

where Y/ m (Fij) are the spherical harmonics correspond- 
ing to the orientation of the vector between particle 
i and its neighbour j in a given coordinate frame. As 
we are interested in local fee-, hep- or rcp-structures, we 
consider I = 6. We assign a vector q 6 (i) to each particle, 
the elements m = — 6 ... 6 of which arc defined as 

06m 00 := -372 ■ W 
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Two neighbouring particles i and j were regarded as 
"bonded", if the dot product q%(i) ■ qh(j) exceeded 0.7 
(i.e. if their local orientational order added up almost co- 
herently). rib(i) is the number of "bonded" neighbours 
of the ith particle. Orientationally ordered clusters are 
defined as regions of bonded particles with nt,(i) exceed- 
ing some common threshold value. In particular, clusters 
were called "crystallites" if rib > 10 (i.e. almost perfectly 
hexagonally ordered). 

Fig. [1] shows system snapshots at two times. Low sym- 
metry clusters (LSC, n& > 5) are light brown, crystallites 
(rib > 10) are green. One can clearly see that the crystal- 
lites are nucleated inside the LSC. To quantify this effect 




FIG. 1: System snapshots at a) 350,000 sweeps and b) 450,000 
sweeps. Particles with nj, > 5 (light brown) and > 10 
(green). Particles with fewer bonds are not shown. 



Fig. [U shows the evolution of crystallinity in more detail. 
In Fig. [5] a) the fractions of particles with given numbers 
of bonds 2 < rib < 12 are plotted. During an induction 
time of ca. 400,000 MC sweeps the amount of orientation- 
ally ordered material grows slowly. Then crystallization 
sets in, as shown by the evolution of particles located in 
pure fec/hep crystallites (rib = 12). The growth rates of 
all other ordered regions are markedly smaller, in partic- 
ular for rib < 5- Fig.[2]b), c) and d) show the evolution of 
crystallites (rib > 10) and LSC (rib > 5) where we distin- 
guished between LSC that contain crystallites and LSC 
that do not. The average number of particles N in (or 
likewise the radius of gyration R g of) the LSC is consis- 
tently larger than N or R g of the crystallites, supporting 
the picture of crystallites growing inside LSC. Further 
analysis showed that the centers of mass of the crystal- 
lites and their surrounding LSC practically coincided. 

To compare the simulation with light scattering data 
we extracted the radially averaged pair correlation func- 
tion as function of time. By Fourier transformation the 
structure factor of the system was calculated. Following 
the procedure first proposed by Harland and van Megen 
111,] we obtained the time evolution of the crystalline 
structure factor (see Fig. [3]) . The simulation results are 
very similar to the ones obtained in experimental inves- 
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FIG. 2: Evolution of crystallinity. a) Fraction X of parti- 
cles in clusters for various values of n D (number of crystalline 
bonds according to Eq. [T] and text thereafter), b) Number 
of particles in crystallites (clusters with rib > 10, open green 
diamonds) and in low symmetry clusters (n > 5) either con- 
taining crystallites (black squares) or not containing crystal- 
lites (red circles), c) Radius of gyration of clusters (symbols 
as above), d) Radius of gyration versus number of particles 
in clusters. 



tigations [20j: At early times during the induction stage 
we observe only one broad peak close to the position of 
the fee (111) peak stemming from the compressed precur- 
sor structures growing slowly in intensity while the width 
remains nearly constant - the scattering pattern changes 
quite slowly. After about 500,000 sweeps a structure fac- 
tor stemming from compressed rhep (random hexagonal 
closed packed) crystals can be identified, the higher order 
reflections are clearly visible: the crystallites evolve from 
the initial precursor structures to a rhep structure. After 
the conversion is complete the structure does not change 
significantly during the main crystallization stage, but 
the intensity increases rapidly (600,000 sweeps till end). 
The peaks become smaller and shift to smaller q values. 

From the integrated area, the width and position of 
individual Bragg peaks the amount of crystallinity, the 
averaged domain size and the volume fraction of the clus- 
ters/crystals can be obtained. While the amount of crys- 
talline material and the lattice constant can be deter- 
mined with high accuracy, there is a noteworthy system- 
atic error in the averaged domain size in our data analy- 
sis, pertaining to Fig- [4j In a polycrystal with rhep struc- 
ture the width of the peaks is connected with the crystal 
size distributions of crystals with different stacking pa- 




q/D 



FIG. 3: Evolution of crystal structure factor. Different col- 
ors correspond to different times given in 10 3 MC sweeps as 
indicated. For further details see text. 
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FIG. 4: Evolution of parameters extracted from analysis of 
the crystal structure factor, a) Crystallinity, red line: data 
from Fig. [2^, (n(,=12) for comparison, b) average domain size, 
c) volume fraction and d) number of clusters/crystals. Black 
squares fee: (111) peak, open circles: fee (220) peak 



rameters [lj|. Especially the analysis of the hcp(002) 
(fcc(lll)) reflection is difficult while the analysis of the 
hep (110) (fcc(220)) peak is quite robust. As we have 
only a very small number of crystals in the sample the 
Scherrer formula to determine the crystal size is strictly 
speaking not fulfilled due to bad statistics. Nevertheless 
we show the determined quantities in Fig. 3] which can 
likely be compared in a qualitative way to experimental 
results [U. 

In the induction stage the amount of orientationally 
ordered material and the size of the highly compressed 
precursors stays nearly constant (2 — 4 x 10 5 sweeps). 
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During conversion the precursors start growing while the 
increase in crystallinity is delayed and a significant drop 
in the number of precursors can be observed (3 — 5 x 10 5 
sweeps). During the main crystallization process also 
information stemming from the fcc(220) peak can be ob- 
tained (5 — 11 x 10 5 sweeps). Here crystallinity increases 
rapidly, crystal size and the number of crystallites show 
their strongest increase due to crystal growth and crystal 
nucleation. The crystals expand to reach the equilibrium 
volume fraction at the end of the crystallization process 
(which was not reached in these calculations). As dis- 
cussed above the absolute values in the averaged crystal 
size stemming from the fcc(lll) are prone to error lead- 
ing to unphysically small values in the absolute number, 
however its time trace reflects the correct trend. 

We would like to emphasize that there is a remark- 
able resemblance of the time evolution of crystallinity 
extracted from the structure factor data with the crys- 
tallinity fraction extracted from the maximally-bonded 
(rib = 12) crystallites (see Fig. 0^). This means that 
the orientation-averaged structure factor analysis can be 
mapped very well to the real-space analysis of the mutual 
particle orientations. 

To summarize, we studied the nucleation process in 
hard spheres by combining a real-space bond-order anal- 
ysis (typical for simulations and confocal microscopy ex- 
periments) with a reciprocal-space analysis of the time- 
evolution of the structure factor (typical for scattering 



experiments). Our simulations showed that nucleation 
in hard spheres is more complex than the traditional pic- 
ture of one-step classical nucleation suggests. We identi- 
fied the formation of dense clusters right after the quench 
containing particles that have high q6q6-coherence with 
at least half of their neighbours. The critical crystal nu- 
clei observed in our simulation are not formed sponta- 
neously in one step from random fluctuations but they 
appear inside these precursors of lower symmetry. The 
metastable fluid relaxes the density first, by producing 
dense low symmetry clusters, and later crystallites of per- 
fect structure are formed. 

The two-step crystallization mechanism identified here 
for hard spheres is akin to processes that have been ob- 
served in protein solutions and in suspensions of attrac- 
tive colloids. Thus we conclude that metastable fluid- 
fluid demixing is not a necessary prerequisite. Different 
dynamics for the two order parameters density and struc- 
ture seem to suffice for a two-step nucleation process. 
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